%% This function estimate the Chen and Szroeter (2012) test
%% Inputs:
% Y=outcome, Z=instrument, D=treatment, set P=1 for step-at-unity P=2 for normal and P=3 for logistic, K=tuning parameter (possible choices sqrt(n/log(n)), sqrt(n/(2*log(log(n))))) for  des=1 if Y is descrete, cl=
% confidence level
%% Outputs
% Rc=1 if the test rejects the null
% Pc=P-value

function [Rc,Pc]=chszpsim(Y,Z,D,P,K,B,A,cl)
%Compute the sample size
n=length(Y);




mu=incosp(Y,Z,D,A);
mu=-mu;
% Estimate the asymptotic variance of sqrt(n)*(mu) (vct is a fuction see the
% related file)

V=vctp(Y,Z,D,B,A);
% Take the main diagonal of the asymptotic variance
theta=1./sqrt(diag(V));

mu=mu(theta<Inf & theta>-Inf & theta ~=0);
V=V(theta<Inf & theta>-Inf & theta ~=0,theta<Inf & theta>-Inf & theta ~=0);
theta=theta(theta<Inf & theta>-Inf & theta ~=0);
switch P
    case 1 % Step-at-unity
        phi=K*theta.*mu<=1;
        % adjustment term
        Lambda =-normpdf(sqrt(n)/K)*(ones(length(mu),1));
        
    case 2 % Normal
        phi=1-normcdf(K*theta.*mu);
        % adjustment term
        Lambda =(-normpdf(K*theta.*mu)*K)/sqrt(n);
    case 3 % Logistic
        phi=(1+exp(K*theta.*mu)).^(-1);
        % adjustment term
        Lambda =((-exp(K*theta.*mu)./(1+exp(K*theta.*mu)).^2))*(K/sqrt(n));
end

% Create a diagonal matrix with the variances of the elements of mu
delta=diag(theta);
% Compute the numerator of the test statistics
Q1=sqrt(n)*phi'*delta*mu-ones(1,length(mu))*Lambda;
% Compute the denominator of the test statistics
Q2=sqrt(phi'*delta*V*delta*phi);
% Compute the P-value
if Q2==0 
    Pc=1;
else
    Pc=normcdf(Q1/Q2);
end
% Reject if Pc<cl
Rc=Pc<cl;